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The distribution and heritability of many traits depends on numerous loci in the genome. In 
general, the astronomical number of possible genotypes makes the system with large numbers of 
loci difficult to describe. Multilocus evolution, however, greatly simplifies in the limit of weak 
selection and frequent recombination. In this limit, populations rapidly reach Quasi-Linkagc 
Equilibrium (QLE) in which the dynamics of the full genotype distribution, including correlations 
between alleles at different loci, can be parameterized by the allele frequencies. This review 
provides a simplified exposition of the concept and mathematics of QLE which is central to the 
statistical description of genotypes in sexual populations. We show how key results of Quantitative 
Genetics such as the generalized Fisher's "Fundamental Theorem", along with Wright's Adaptive 
Landscape, emerge within QLE from the dynamics of the genotype distribution. We then discuss 
under what circumstances QLE is applicable, and what the breakdown of QLE implies for the 
population structure and the dynamics of selection. Understanding of the fundamental aspects of 
multilocus evolution obtained through simplified models may be helpful in providing conceptual 
and computational tools to address the challenges arising in the studies of complex quantitative 
phenotypes of practical interest. 
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I. INTRODUCTION 



R.A. Fisher's celebrated "Fundamental theorem of natural selection", relating the rate of change in the average 
fitness to the genetic variance in fitness, occupies a place in Population Genetics similar to Newton's U F - 



Physics. Yet conceptually Fisher's law and the whole subject of "Quantitative Genetics" (Falconer and Mackay 1996 



Lynch and Walsh 1998), which studies the response of quantitative traits to selection, is closer to Thermodynamics. 
Thermodynamics is a phenomenological description of readily measurable physical properties (e.g. average energy 
or pressure) of a large ensemble of molecules. Quantitative Genetics is a phenomenological description of readily 
observable phenotypic traits of a population. Thermodynamics takes macroscopic averages over the random motion 
of individual molecules in thermal equilibrium. Quantitative Genetics similarly focuses on the behavior of population- 
wide averages (and variances) over many genetically diverse individuals. The genetic composition of the population 
is governed by natural selection and random drift along with recombination and mutation, all acting on individuals. 
The phenotype distribution is related to the genotype distribution by the largely unknown genotype-to-phenotype 
map, which is further obscured by environmental effects which can cause phenotypic variation even between genet- 
ically identical individuals. Yet deterministic laws of Thermodynamics emerge despite the complexity and chaos of 
molecular motion. In fact they emerge thanks to the microscopic complexity and chaos and are made possible by the 
extensive self-averaging that dominates macroscopic behavior of physical matter. Similarly, simple laws of quantitative 
population genetics emerge when phenotypic traits depend on large numbers of polymorphic genetic loci. 

While the analogy between Quantitative Genetics (QG) and Thermodynamics is most appealing and has been noted 

see 



1930) 



Barton and Vladar (2009); Iwasa (1988); Sella and Hirsh 



by many including R.A. Fisher himself (Fisher 
(2005) for recent work - fundamental issues such as the lack of energy-like conserved quantity in population genetics 
impede direct transcription of thermodynamic laws to QG. Instead, the analogy must be pursued as an approach 
to the construction of a coarse-grained phenomenological theory bridging the gap between ensemble averaged (read 
population averaged) observables and the hidden micro-scale (read individual genotype) dynamics. One must be 
careful to define an averaging ensemble that equilibrates on the time-scale of the observation, e.g. the response 
to selection in QG. In particular, as we illustrate in Figure [I] dynamics in sexually reproducing populations are 
characterized by two widely different time-scales: 1) mating and recombination reshuffle the polymorphic loci, allowing 
exploration of the space of genotypes on a short time scale and 2) mutation and population drift control genetic 
variation on much longer time scales, often long enough to render the ensemble meaningless. 

The bridge between the dynamics of the genotype distribution and the coarse grained, QG-type, description is 
built on understanding multi-locus evolution. Our review will focus on the intermediate time scale in the above 
mentioned hierarchy. We will show how the genotype distribution P(g,t) can be parameterized by slowly varying 
allele frequencies, while mating and recombination lead to rapid equilibration of P(g, t) given a set of allele frequencies. 
In this ensemble trait distributions are determined by allele frequencies and the dynamics of trait averages can be 
expressed in terms of the dynamics of allele frequencies. This in turn gives rise to the familiar laws of quantitative 
genetics in terms of additive variances and covariances. In this sense, a statistical multi-locus theory plays the role of 
Statistical Mechanics, which explains how the deterministic laws of thermodynamics emerge from the erratic motion 
of many microscopic particles. Hence the subject of the present review should be thought of as "Statistical Genetics" 
- a term introduced in a closely re lated context by Wright (1942). 

Classical Quantitative Genetics (Falconer and Mackay 1996) is based on the assumption that genotypes are random 
re-assortments of alleles, each occurring with a certain frequency. This absence of correlations between alleles at 
different loci is termed "linkage equilibrium" , implying that recombination (breaking linkage) has relaxed correlations 
between loci. This drastic simplification has earned QG a derogatory epithet of "beanbag genetics" from the pen of 
Ernst Mayr (MayrJ 1963) (see however Haldanc (1964) in defense of beanbag genetics). Yet in the present review 
we shall see that the key phenomenological laws of QG extend beyond the assumption of linkage equilibrium. This 



understanding emerges from the studies of multilocus selection which began with two alleles/two loci systems (Karlin 



and Feldman, 1969 Kimura 1956| Lewontin and Kojima 1960). 



Kimura ( 1965 1 showed that a two locus system 
tends towards a state where allele frequencies change slowly and correlations are small and steady. He termed 
this state Quasi-Linkage Equilibrium (QLE), which is the subject of this review. Subsequently, several comprehensive 
treatments of multilocus evolution were developed ( Baake 2001 Barton and Turelli 1991 Burger 1991 Christiansen 



Nagylaki 1993| |Prugel-Bennett and Shapiro 1994 j(for a monograph see|Burger ( 2000 )) with Barton and Turelli 



1990 



( 1991 ) and NagylakiJ ( |1993[ ), in particular, generalizing and justifying the QLE approximation in multi-locus systems. 



In addition to the study of generic behavior of systems with a very large number of loci, explicit multi- locus modeling 



of smaller systems has been used to study the evolution of recombination (Barton 



patterns o f genetic variation produced by positive selec tion (|Stephan et aL[|2006|) . Recent work produced interesting 

fitness landscapes with five or more 



2009 Weinreich et al. 



2006 ) of empirically determinec 



1995a Roze and Barton 2006 1 and 



examples ( de Visser et al. 

loci. The dynamics of populations on these landscapes can be studied in laboratory experiments and comparison to 
theoretical models is possible (de Visser et al. 2009). Quantitative understanding of multi- locus evolution is also 



3 



Long times : Ensemble of Populations 
e.g. distributions of allele frequencies 

• Kimura's diffusion equation 

• Mutation, selection, drift balance: 
Wright's Equilibrium 




t 



Mutation, drift, selection over long times 



Intermediate times: Ensemble of individuals 


1 








• Allele frequency dynamics and 


^2 


_ — y 1 
, _ ~ y / 1 




Quantitative Genetics 


, - , s 1 1 




• Wright's adaptive surface 




. . « i 1 




• Fisher's Theorem 










u o 


1 



t 



Recombination and 
selection on standing variation 



Short times : Dynamics of genotype distribution 

• "Microscopic" description: 

Pfat), 9 = {*!,-. -,sl} 

• P(g,t) is defined on 2 L possible genotypes 



- + 



+ 



1® 



FIG. 1 Time scales in sexual populations. A population is described by the distribution of 2 genotypes, on which selection, 
recombination and mutation acts. In sexual populations, mating and recombination is the fastest process, so that different loci 
are only weakly correlated (close to linkage equilibrium) and its dynamics can be approximately described via L allele frequencies 
- a number much smaller than the 2 L possible genotypes that would have to be tracked otherwise. Allele frequencies change 
slowly and means of quantitative traits follow the laws of quantitative genetics. Over the much (much) longer time-scale of 
/x -1 , allele frequencies themselves tend to an equilibrium between selection, mutation, and genetic drift (assuming a constant 
environment). 



essential when studying the emergence of drug resistance in HIV, which often depends on several interacting loci in 



a recombining population ( Bretscher et al. 2004 



Our discussion of the multilocus selection pro 



Kellam et al. 1994 Nora et al. 2007). 



jlem will follow the Barton- Turelli course of making and keeping 
it simple (Barton and Turelli 1991 Turelli and Barton 1994), by attempting to make it simpler still. We will 
define a streamlined conceptual and analytic framework which will not only reproduce classic results, but also readily 
generate some new extensions. To accomplish this we will formulate and analyze a "minimal model" of multilocus 
evolution: continuous time selection in a haploid model with a general (epistatic) fitness function of L loci. We shall 
show how the Generalized Fisher Theorem and other results of quantitative genetics follow from a straightforward 
cumulant perturbation theory similar to that used extensively (for high temperature expansions) in Statistical Physics 
(McQuarrie 1973). The perturbative regime corresponds to Kimura's Quasi-Linkage Equilibrium (QLE). Using this 
formulation of QLE, we present systematic generalizations of QG results and of (Kimura's) diffusion theory, typically 
formulated in complete linkage equilibrium, to include weak correlations between loci. We also discuss how QLE 
breaks down when the ratio of characteristic strength of selection to the rate of recombination exceeds a critical 
value that depends on the strength of epistasis. While the QLE regime corresponds to selection of individual alleles 
based on their effect on fitness averaged over genetic backgrounds, the breakdown of QLE follows the appearance of 
strong correlations between alleles at different loci and represents a transition to effective selection of genotypes. In 
the Discussion section we shall connect the transition from "allele selection" to "genotype selection" to the closely 



related spin-glass transition (modeling the behavior of disordered magnets) studied in Statistical Physics (Mezard 



et al. 1987). We shall also discuss its implications for Quantitative Genetics. 
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II. RELATING QUANTITATIVE TRAITS AND GENOTYPES. 



Let us focus on the fitness which is the most important example of "quantitative trait" , although everything we 
shall say about "fitness landscapes" in this section, applies directly to any quantitative phcnotype. A fitness landscape 
is a metaphor for a map from the high dimensional space of genotypes to expected reproductive success. While the 
map itself is unambiguous, several different ways of parameterizing fitness landscapes with alleles and groups of alleles 
have been proposed (Barton and Turelli 1991 Hansen 2006| Hansen and Wagner 2001 Weinberger 19911. 



Consider a haploid genome of L loci with two alleles each, such that a genotype is uniquely characterized by L 
binary variables g = {si, . . . , sj,}. We choose G {—1,1}, i = 1,...,L instead of Sj £ {0,1} (more commonly 
used in population genetics literature), since the symmetric choice simplifies the algebra below. (The relation between 
representations can be found in Appendix [5] a short glossary of population genetics terminology is given in Appendix 
|A| ) Functions of the genotype, e.g. as the population distribution, fitness, or any other quantitative trait, live therefore 
on a L-dimensional hypercube. Any such function on the hypercube can be decomposed into a sum of monomials in 



S i S j + E fijkSiSjSk + 
i<j<k 



(1) 



where the first sum represents independent contribution of L single loci, the second sum which runs over all L(L — 1) /2 
pairs of loci represents contribution of pairs and the higher order terms account for the effect of each and every possible 
subgroups of loci. The first order contribution fi defines the additive effect of locus i which is independent of all other 
loci considered. Higher order terms which include locus i define the genetic background dependence of the effect of 
Si allele. Collectively, terms of order higher than one represent genetic interactions also known as "epistasis". The 
contribution of each locus or subgroup of loci is determined by unbiased (i.e. each genotype enters with weight 2~ L ) 
averaging over the remainder of the genome: thus the coefficients are given by 



(2) 



One easily convinces oneself that plugging Eq. ([lj into the expressions in Eq. ^ reduces to the desired coefficients. 



In total, there are 2 L coefficients ik , as it has to be for an exact representation of a function on a hypercube. 
In fact, the coefficient of the expansion of F(g) into monomials is nothing but the Fourier transform of the original 
function on the hypercube, which was used in the context of genotype-fitness maps in Hordijk et al. (19981; Stadler 



and Wagner (1998); Weinberger (1991). In addition to this genetic contribution to the trait, the trait value of a given 



individual will also depend on environmental (and epigenetic) factors which are not modeled here. 

It proves very useful to define a "density of states" p(F) = 2" L ][] g <5(F - F(g)), where 8(F) is a Dirac delta- 
function. The fraction of genotypes with fitness in the interval [F, F + 8F] is then given 

byjf dF'p(F'). Provided 

F(g) receives contributions of very many terms of similar magnitude in Eq. [TJ the Central Limit Theorem will apply 
making the density of states approximately Gaussian in shape. The width of this Gaussian is given by the (square 
root of the) variance of F (g) over the hypercube: 



= 2 



-L 



1.1 k 



(3) 



i<j 



i<j<k 



This simple decomposition of variance is the equivalent of the Parseval's theorem for the Fourier transform. Note 
that this variance is an intrinsic property of the fitness landscape completely independent of any population that may 
be evolving on it. It should not be confused with the population variance that we will discuss later. We shall use a 
as a measure of selection strength. 

The sums in Eq. (§ for a can be interpreted as the power spectrum of the F(g). A falling or rising power spectrum 
gives rise to qualitatively different landscapes: If most of the variation of the fitness function were captured by the 
first order terms, the landscape would be smooth and simple. If higher order terms dominate the fitness variance, 
the landscape is multi-peaked and rugged. The properties of smooth versus rugged landscapes (parameterized in the 
manner of Eq. ( fl|) ) h as been a subject of extensive study in statistical physics as it relates to the theory of spin-glasses 
(Mezard et al. |1987). It is known that the key consequences of complexity of the general landscape, appear already 



in the class of functions involving only pairwise interactions (a.k.a. the Sherrington-Kirkpatrick model (Sherrington 
and Kirkpatrick 1975)) (Mezard et al. 19871. Here for simplicity we shall consider only pairwise interactions. (An 
alternative instructive simplification would be to consider F (g) defined by a fixed random function on the hypercube, 
known in population genetics as the "house-of-cards" model (|Kingman 1978 ) or NK-models Kauffman and Weinberger 
(1989) and in physics as a "random energy" model ( Derrida[ " 1981 \7) 
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Before moving on to population dynamics, it is instructive to discuss the implication of the combinatorial explosion 
of higher order interactions: In principle there are (^) interactions of order k, a number which increases with L as 
L k . Hence increasing the number of loci without changing the statistics of the coefficients would shift the power 
spectrum towards higher order, making the function more rugged. It seems more likely that the interactions are 
sparse with the number of "partners" of a typical locus not growing in proportion to the total number of loci: In 
particular one may posit that each locus interacts with a finite number of other loci, independent of L and set all 



other coefficients to zero. Unfortunately, despite some recent progress flBrem and Kruglyak 2005 



Ehrenrcich et al. 



2010), we still know very little about generic structure of genotype-phenotype maps. One must also be aware of the 
fact that, because of selection, statistics of genetic interactions observed among co-segregating polymorphisms within 
a breeding population may be quite different from that for a random set of loci or a for polymorphisms created by 

Indeed, most immediate evidence for epistasis is provided by the 



crossing two isolated populations ( Jinks et al 
"outcrossing depression" 



1966). 



1966 Seidel et al. 2008). 



suppression in the fitness of progeny issuing from a cross of diverged strains ( Jinks et al. 



III. DYNAMICS OF THE GENOTYPE DISTRIBUTION. 

Selection, mutation, and recombination operate on individuals and change the distribution of genotypes, P(g, i), in 
the population. The fitness F(g) of a genotype g is defined as the expected reproductive success, i.e. the rate at which 
the proportion of a genotype increases or shrinks in the population due to (natural or artificial) selection. During the 
time interval At, selection changes the distribution of genotypes according to 

AtF(g) 

P(g,t + At) = -^ MF jP(g,t) (4) 

where (e AtF ) — ^2 g e AtF ^ P(g,t) denotes the population average. The genetic diversity that selection acts upon is 
due to mutations, which change the genotype distribution as follows 

L 

P(g, t + At) = P(g, t) + Atn £ [P(M i9 , t) - P(g, t)] . (5) 

i=l 

Here M^g is a shorthand for genotype g with Sj replaced by — Sj. Despite the importance of mutations for generating 
polymorphisms and maintaining genetic diversity in the long run, the effect of mutation on the dynamics of significantly 
polymorphic sites can be neglected if mutation rates are much smaller than selection coefficients. 

In addition to selection and mutation, the dynamics of the genotype distribution in sexual populations are driven 
by mating and recombination. Gametes are formed during meiosis crossing over homologous parental chromosomes. 
Assuming random pairing of gametes and outcrossing with rate r, the genotype distribution changes during recombi- 
nation as follows: 

P(g,t + At) = (l-Atr)P(g,t)+Atr £ C({Q)P(gW ,t)P(g^ ,t) . (6) 

{&}«} 

The first term accounts for those individuals that did not outcross during the At time interval. In the event of 
outcrossing, a new genotype is formed from genetic material of the mother with genotype g^ 1 and a father with 
genotype g^\ The novel recombinant genotype g inherits a subset of his loci from the mother and the complement 
from the father, which in Eq. ^ is described by the set of random variables {^}. If & = 1, gene i is inherited from 

the mother, if £j = from the father. Using this notation, the maternal genotype is s^™ 1 ' = £jSj + (1 — £i)s^ and 

( f) 

equivalently the paternal genotype s { = ( 1 — £j ) Sj + & s 'i ■ The part of the maternal and paternal genome which is not 
passed on to the offspring, {s^}, is summed over. Each particular realizations of i.e. a pattern of crossovers, has 
probability C({£}), which depends on the crossover rates between different loci. In addition to the summation over all 



{s-}, we have to sum over possible crossover patterns {&}. A very similar notation was used in Christiansen (19901. 
While our presentation so far was completely general, dealing with diploid genomes inflates the required book-keeping 
as we proceed with the analysis. Since our goal is to present the key effects and ideas in the simplest possible form, 
we shall from here on restrict to considering only haploids, two of which recombine upon mating producing a haploid 
offspring. Although this model is chosen for simplicity sake, it is sufficient to describe diploids in the absence of 
dominance. It also describes haploid yeast going through mating/sporulation/germination cycle or the population 
genetics of many RNA viruses like HIV and Influenza. 
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Provided selection is weak (AtF(g) <C 1), we can use a continuous time description of the dynamics. 

l 



dt 



P(g,t) =(F(g) - (F))P(g,t) + (xJ2l P ( M i9) ~ P($)] 

i=l 

+ r £ C{{$)[p{g^\t)P{gU\t)-P{g',t)P{g,t) 



(7) 



This equation describes the dynamics of the genotype distribution in the limit N —> oo where each genotype is sampled 
by enough individuals to neglect sampling noise which would arise during reproduction. This stochastic component 
to the dynamics of the genotype distribution is known as "random genetic drift" . We shall discuss random drift 
in Section |VI| Our focus here will be on the interplay between selection and recombination, which dominates the 
behavior of Eq. 0. 

Instead of specifying P(g, t) for every g, P(g, t) can be parameterized by its cumulants. The cumulants of first and 
second order are defined as Xi = (sj) and Xij = ( s i s j) ~ ( s i)( s j)i which are related to allele frequencies and pairwise 
linkage disequilibria (see Table Ilk. In total there are 2 L — 1 cumulants, with higher order ones, Xij...ki more easily 
defined via the cumulant generating function ( McQuarrie| |1973| . However, only the first and second order cumulants 
will be needed in the present context. 

To obtain dynamical equations for Xit we multiply Eq. |7| by Sj and sum over all possible genotypes. One finds 



Xi = (s i [F(g)-(F)])-2fx(s i ) , 



(8) 



where we have used MjSj = —Si and used the notation Xi f° r total derivative with respect to time. The dynamics 
of Xi do not depend explicitly on the recombination rate, which is intuitive since recombination does not create or 
destroy alleles. In order to evaluate (siF(g)) in Eq. ^ we need to know higher order cumulants, i.e. we are faced 
with a hierarchy of cumulant equations. 

In contrast to first order cumulants, the dynamics of higher order cumulants depend explicitly on recombination, 
which has the tendency to destroy associations between alleles and drives higher order cumulants to zero. To write 
down an equation for the dynamics of the second order cumulants, Xij — ( s i s j) ~ XiXj 
which explicitly depends on recombination. Evaluating the recombination term only, we find 



we have to evaluate -^(siSj), 



■£<?({£}) 



P(g^\t)P(g^,t)-P(g',t)P(g,t) 



CijXij 



(9) 



where cy is the probability that loci i and j derive from different parents: Cjj = — £j) + (1 — 

To arrive at this result, we substituted Sj = £is| m ' > + (1 — (analogously for Sj), and averaged over the maternal 

and paternal genomes. The second term evaluates simply to r(siSj). This result holds more generally for central 
moments of the genotype distribution ( Barton and Turelh] 1991[ ). Together with selection and mutation, we find (for 



Xij = ((si ~ Xi)(sj ~ Xj)(F(g) - (F))) - 4nxi 3 - rc^j 



(10) 



We see that selection drives Xij away from zero, while Xij relaxes though mutation and recombination. In absence of 
selection P(g, t) tends to a steady state of "linkage equilibrium" (LE) with vanishing cumulants Xij (f° r * 7^ j) implying 
complete decorrelation of alleles at different loci corresponding to factorization of the genotype distribution: Pp(g) = 
Ylf—i Pi(si)- It is easy to see that the recombination term in Eq. |7| vanishes whenever P(g) — Po{g)- In Section|v[ we 
will, starting from Po{g), develop the Quasi-Linkage-Equilibrium (QLE) approximation by systematically accounting 
for small linkage disequilibria (Xij)- 



IV. TRAIT DISTRIBUTIONS AND THE DYNAMICS OF POPULATION AVERAGES 

In most cases, P(g, t) cannot be observed directly. Instead, the subject of quantitative genetics are distributions 
of traits in the population. Trait distributions can be obtained from genotype distributions by projection. The 
probability of finding in the population an individual with fitness (or any other trait) in the interval [F, F + AF] is 
given by 



p(F,t) = ^28(F-F(g))P(g,t) 
a 



(11) 
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where 6(F) is the Dirac delta-function (J dF6(F) = 1). Applying this projection to Eq. ^ yields an equation for the 
dynamics of the trait distribution. Before addressing the dynamics of traits in sexual populations, it is instructive to 
consider the dynamics of the fitness distribution p(F, t) in absence of mutation and recombination, in which case one 
obtains simply 

j t p(F,t) = [F-(F)]p(F,t) (12) 

where (F) = J dFFp(F,t). Multiplying this equation by F and integrating over F (i.e. the 1st moment of this 
equation) yields Fisher's "Fundamental Theorem" in the asexual case. 



dt 



(F) = ([F (F)} 2 ) 



(13) 



Evidently this is just the 1st i n th e hierarchy of infinitely many moment equations that characterize the dynamics of 
p(F,t) given explicitly by Eq. (12 1. The 2nd moment expresses the dynamics of a 2 in terms of the 3rd moment, etc. 



This hierarchy of equations is not closed, yet under certain co ndi tions higher moments may be suppressed making 
the a 2 a slowly varying function of time. One notes that Eq. ( 12 ) has a Gaussian traveling wave solution p(F, t) = 
Cexp[— (F — vt) 2 /2v] with an arbitrary constant variance a 2 = v setting the rate of fitness growth j^(F) — v in 
agreement with Eq. ( |13[ ) . A traveling wave with constant speed requires that genotypes with arbitrarily high fitness are 
populated with at least one individual, which requires an infinitely large population with infinitely many polymorphic 
loci with limits taken in this order. Otherwise genetic diversity disappears and adaptation stalls. The evolution of 
the shape of the fitness distribution in finite populations has been studied in the context of Genetic Algorithms by 
Prugel-Bennett and Shapiro ( |1994 ) . Priigel-Bennett and Shapiro study the effect of selection and recombination on 
the cumulants of the fitness distribution and observe how the fitness variation vanishes as the population condenses 
into a local fitness maximum. To prevent this condensation, new variation has to be constantly supplied by mutation. 
Quite generally a is determined by the balance generation of genetic variation through mutations or recombination and 



its removal be selection and drift, which requires careful stochastic treatment (Desai and Fisher 2007 Hallatschek 



2011 Neher et al. 2010 Rouzine and Coffin 2005 Rouzine et al. 2003 Tsimring et al. 1996 1 



One can also consider the dynamics of an arbitrary trait G(g) different from fitness. In analogy to Eq. (11 1, we can 



study the joint distribution, p(F, G, t) of this trait with fitness. The population average of the trait obeys 

d 



dt 



(G) = (GF) - (G)(F) = Cov(F,G)p( g ) 



(14) 



i.e. its rate of change is given by its covariance with fitness ( Price| 1970[ ). This statement is also known as the 



"secondary theorem" of natural selection (Robertson 1966) 



With mutation and recombination, the dynamics of trait means are no longer that simple. To evaluate the mutation 
and recombination terms, we utilize the orthogonal expansion of the fitness function in Eq. ([!]). Restricting ourselves 
to pairwise interactions, we can use Eqs. ^ and Q to obtain 

— (F) = a 2 - fiAf, -r^2 c ijfijXij > (15) 

where A M is the average loss in fitness due to mutation. The latter can be calculated by observing that each moment 
decays through mutation with rate 2/iA:, where k is the order of the moment. 



2 /** + 4 bbXj + Xij ) 



(16) 



Higher moments decay faster because they have a greater mutation target. The second term in Eq. (15 1 is the loss 



in fitness through recombination, which reflects the tendency of recombination to factorize the genotype distribution 
such that contributions like fijXij to (F) decay with rate rcij. We will later see that the previous form of Fisher's 
theorem can be recovered by a suitable definition of an additive fitness variance. To do so, however, we have to 
understand how the genotype distribution evolves under selection and recombination. 



V. BEYOND LINKAGE EQUILIBRIUM: QUASI LINKAGE EQUILIBRIUM 

We have already seen that without selection or without epistasis, P(g) distribution will relax to a product of 
independent distributions at different loci: the linkage equilibrium state. Next we shall account for the correlations 
between loci induced by selection. For simplicity we shall omit the mutational contribution, which we shall restore 
once we understand the basis of Quasi Linkage Equilibrium (QLE). 



A. QLE: A perturbation expansion at high recombination rates 



If selection on the time scale of recombination is weak, i.e. 



be calculated using perturbation theory (Barton and Turelli 
genotype distribution as follows 



a <C r, the induced correlation is also weak and can 
1991 Kimura 1965). To this end, we parameterize the 



log P{g,t) =$(*)+ Mt>i 



i<j 



(17) 



which is the already familiar Fourier representation of functions on the genotype space. The factorized distribution 
Po(g) would correspond to the coefficients, 4nj, of all multilocus contributions being zero. The second order terms 
capture (to the leading order) the correlations induced by selection and (in the limit under consideration) are assumed 
to be small. The genotype independent term $(£) is fixed by the normalization of the probability distribution. 



-*({«) _ 



E 



cxp 



<piSi + (f>ijs l s j 



(18) 



and acts as the generator of the cumulants via 



<9$ 



<9 2 $ 



dcj)i ' lJ d<f>id(f>j 

The generating function $ is evaluated perturbatively for small <pij in the Appendix [C] yielding 

Xi 



tanh(0i) + 4>ij0- ~ tanh 2 (</>i)) tanh(^j) 



Xij 

Xu 



l-Xi 



for i ^ j 



(19) 



(20) 

(21) 
(22) 



which is correct to the leading order in \<f>ij\. The distribution given by Eq. (17 1 may be thought of as a maximum 
entropy distribution constrained to have certain first and second order cumulants: Parameters (pi and ipij are the 
Lagrange multipliers that impose the constraints. 

Let us rewrite Eq. ^ as an equation for the dynamics of logP(<?) which yields 



P{g {m) )P(g {f) ) 



P(g)P(g / ) 



i<3 



- 1 



(23) 



where the recombination part has been evaluated approximately by expanding the exponential that defines P{g) (see 
Appendix [C]) . We can now collect terms with the same monomials in Si to obtain the equations governing the time 
evolution of d>j and <&,•: 



— fij TCijtfiij 



(24) 
(25) 

At large crossover rates rc^, the (fry rapidly approach a steady-state <pij — i^r.- This has to be contrasted with the 
behavior in absence of recombination, in which case <f>ij would grow linearly as fyt. Recombination prevents effective 
selection on interactions. Instead, the higher order contributions to fitness affect the dynamics of after averaging 
over possible genetic backgrounds: Substituting the steady state relation into the equation for yields 



<j>i = fi + E -fijX] = fi 



(26) 
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where we have defined fi = fi + J2j fijXj which is the effective strength of selection acting on locus i in linkage 
equilibrium. It is obtained from the general expression for F(g) in Eq. ([I]) by replacing Sj ~ > Xj an d differentiating 
with respect to Xi (and is truncated here at second order because we assumed, for simplicity, that genetic interactions 
are limited to that order). 

Converting <ps to X s using the relation (20), we find Xi — (1 ~Xi) ( Pii correct to the leading order. For the discussion 



below, it will be useful to derive equations for Xi an d Xij a ^ so to the sub-leading order 



3 



Xij 



fj-Xifij +<tO(a 2 /r 2 ) 



(27) 



Xij 



(i-xf)(i-x J 2 )fe 

ViX.% + 2fjXj + rcij 



+ 0{a 2 /r 2 ) for i ^ j 



In QLE, correlations Xij between loci (i ^ j) are determined by the balance between epistatic selection and re- 



combination, 
frequencies.) 



(Note, in contrast, the diagonal elements Xi 



(s 2 



[Si)' 



1 — xt are determined by the allele 



Wright ( 1931 ) showed that in linkage equilibrium, the dynamics of allele frequencies are driven by the gradient in 



mean fitness. The result can be generalized to include correlations between loci arising in QLE. Starting with the 
exact equation for the allele frequency dynamics and using our parameterization of P(g) via the "fields" <f>i given in 
Eq. JPTl), we find 



Xl = (a i F)- X i(F) = d^(F) « J2 d ^d Xj (F) 

j 



(28) 



where we have used the chain rule of differentiation and the fact that = Xij following directly from Eq. (19). The 
correlation matrix Xij ac ts as a mobility matrix for allele frequencies. The non-diagonal entries of order a jr imply 
that selection on locus j, via the correlation with locus i, affects the rate of change of Xi- Eq. (28) describes the 
dynamics of allele frequencies as the population ascends Wright's "adaptive landscape" . While allele frequencies still 
evolve to maximize (F), their dynamics now are coupled by correlations captured in the off-diagonal terms of Xij- 

The key point emerging from the analysis of the weak selection/rapid recombination limit is the remarkable simplic- 
ity of multi-locus dynamics: the 2 L ordinary differential equations for all cumulants or equivalently for all genotypes 
are reduced to L differential equations describing the dynamics of allele frequencies. Higher order cumulants are 
slaved to allele frequencies and can be obtained by solving algebraic equations defining the L dimensional quasi- 
linkage-equilibrium manifold. The distribution of genotypes in the population can therefore be parameterized by 
time-dependent allele frequencies, with all other features of the distribution constrained by the QLE equations. In 
mathematical terms, the dynamics of genotype distribution are approximately reducible to the dynamics on the "center 
manifold" formed by the set of allele frequencies ( Guckenheimer and Holmes 1997 ) . Within the QLE approximation, 
population averages of any trait G(g) can be parameterized by {xi{t)i ■ • • j XiJtjJ an d the time-derivative of the trait 
mean is therefore given by 



dt 



(G) ^a x ,(G)5 fX ,(t) =Y.xMG)d Xj (F) - 2^ Xi S Xi (G) 



(29) 



where we have restored the contribution of mutations through its effect on allele frequencies as it appeared in Eq. {[sj) . 
This result has a very simple interpretation: The rate of change of the trait mean is the product of the rate of change 
of allele frequencies through selection and the susceptibility of the trait mean to the allele frequency. The second term 
accounts for the effect of mutation on the trait mean. Since the first term is the additive covariance between fitness 
and the trait G, this equation is the analog of Eq. ( 14 ) in a recombining population. The QLE approximation breaks 



down when recombination is not sufficiently rapid to confine the genotype distribution to the L dimensional manifold 
defined by quasi-steady correlations between loci. This breakdown will be discussed in more detail below. 



B. Additive genetic variance and Fisher's theorem in QLE 



Fisher's theorem in sexual populations posits that the rate of mean fitness increase is equal to the additive variance. 
We will now discuss how Fisher's theorem emerges from Eq. ( 15 1 and how it compares with Eq. (29 1 which obviously 
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can be used to calculate d(F) /dt. Additive variance is typically denned as the variance captured by a linear model of 
the form 



^4(3) = ao + ^aiSi 



(30) 



where the coefficients are determined by minimizing 

a} = Y,{FA{g)-F(g)) 2 P{g,t) 



(31) 



The remaining variance o~ 2 is commonly called epistatic or interaction variance. Minimization yields ao = (F)—'^2 i a,i\i 
with <2j determined by the linear equation 



Y,Xi j a j = {s i F)-Xi{F)=d 4>i {F) 



(32) 



We have seen the right hand side of this equation already in Eq. ( 28 ) : it is the contribution of selection to \i . In the 
high recombination limit, d ( j >i (F) « J^j Xij^xj (F)- Hence the additive fitness coefficients (defined by linear regression) 
are <2j = d Xi (F), which is accurate to order a jr. The additive variance therefore is 



a iXij a j 

ij 



Y,xA(F)d X] (F)+a 2 0(a 2 /r 2 ) 



(33) 



Recalling the QLE equation for mean trait dynamics, Eq. (29), and using fitness as a trait, we have 

d 



dt 



(F) « £ Xijd Xi (F)d Xj (F) -2^ Xid Xi (F) 



(34) 



and comparing to the definition of a\ we arrive at the generalized Fisher's "Fundamental Theorem" 



dt 



(F) 
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^ + 0(a 4 /r 2 



(35) 



which limits growth of fitness to the additive variance. Comparing to the general expression for mean fitness given 
before in Eq. ( 15 ) we see that the loss in fitness due to disruption of favorable combinations of alleles though 
recombination exactly cancels the epistatic a 2 = a 2 — a\ part of total variance. In other words, in a sexually 
reproducing species the uncertainty in the phenotype of the offspring in relation to that of its parents limits the effect 
of selection to the additive component of variance. The latter is that genetic component of the trait that "survives" 
reshuffling of genes by rcassortmcnt and recombination which depends on the genetic distance to the mate. Hence, 
this decomposition of genetic variation in additive and non- additive components is explicitly population dependent. 

One must of course remember that the generalized Fisher's law as stated only holds in this rapid recombination/ weak 
selection limit and only after correlations have relaxed to their steady QLE values. During the initial transient towards 
QLE or at low recombination rates mean fitness can exhibit very different dynamics. The meaning of Fisher's theorem 



has been subject to extensive discussion in the literature (Edwards 1994 Ewens 1989 Feldman and Crow 1970 



Frank and Slatkin |1992 Price 1972) caused by Fisher's insistence that his statement was exact. Price " p~972 ) in 



particular suggested that Fisher's intention was to describe not the total rate of change of mean fitness, but only 
the "partial rate" due to change in allele frequencies: i.e. just 1st term on the r.h.s. of Eq. (29). The "theorem" 
would in that case become an exact statement, but not a very useful one! Following Kimura ( 1958 ) and Nagylaki 



(1993) our Eq. (35) sticks to d(F)/dt so that the generalized Fisher's theorem is an unambiguous, but approximate 
statement. The above analysis assumed that the population is subject to a constant fitness function and the mean 
fitness provides a useful measure of adaptation. If the fitness function itself depends on time, the increase in mean 
fitness due to adaptation of the population is superimposed with the dynamics of the fitness function. In the latter 
case, an unambiguous measure of adaptation, the fitness flux, can be defined in analogy to fluctuation theorems of 
non-equilibrium statistical mechanics ( |Mustonen and Lassig| ["2010[ ). 

The off-diagonal terms in the additive variance aiXijdj have interesting implications for the evolution of recom- 
bination: if two alleles that are selected with the same sign (e^a., > 0) are anti-correlated < 0), the rate of 
adaptation is smaller than it would be in linkage equilibrium. This is the basis for the often made statement that 
recombination accelerates adaptation by reducing negative linkage disequilibria and thereby increasing the additive 
variance (Barton and Otto 2005). There is, however, an additional effect of recombination on adaptation that is not 
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captured by deterministic multilocus dynamics and is likely to be more important: Recombination greatly increases 



the likelihood that a novel beneficial mutation establishes and ultimately fixates in the population (Barton 1995b 



Fisher 


1930 


Muller 


1932 


Neher et al. 


2010) 



Thereby the number of simultaneously polymorphic loci is increased, 
which in turn increases the fitness variance and speeds up adaptation. The reason for this is again that recombination 
breaks down negative linkage disequilibria (a tendency of beneficial alleles to be anti-correlated) , which are generated 

Analysis of this phenomenon requires going beyond 



by chance and amplified by selection (Barton and Otto 
QLE (see below). 



2005). 



VI. FINITE POPULATION DRIFT AND WRIGHT'S MUTATION/SELECTION/DRIFT EQUILIBRIUM. 



So far our formulation of the genotype dynamics Eq. Q and the dynamics of allele frequencies Eq. p7| was 
deterministic, i.e. we neglected random drift. Random drift is a consequence of the stochastic nature of birth and 
death in a finite population of size N . In the simplest models of stochastic population genetics - called Fisher- Wright 
models - stochasticity is introduced by resampling the population from a multinomial distribution parameterized with 
the current genotype (or gamete) frequencies each generation. 

We have seen above that the genotype frequency distribution can be parameterized by allele frequencies when 
recombination is rapid and we shall discuss now how resampling of genotypes leads to stochastic contributions to 
the dynamics of allele frequencies and cumulants. For alleles that are present in large numbers, the relative sizes of 
fluctuations due to resampling are small and random drift can be accurately described by a diffusion approximation 
( Kimura 1964 ) . To derive a diffusion equation for allele frequencies, we generalize the ordinary differential equations 
Eq. (27) to stochastic differential equations (Langevin equations (Gardiner, 2004 1). For a finite time step At, one has 



Xl (t + At) = x<(t) + At 



^AtQ(t) 



Xij (t + At) = Xij {t) + At [(1 - - X 2 j)fij - r^] + VAtQj(t) 



(36) 



(37) 



where we have neglected terms much smaller than rc^ in the relaxation rate of Xij ■ G (t) an d Cij (t) are white noise 
terms with zero mean and a covariance matrix determined by the multinomial sampling of the genotypes. One finds 



(Ci(t)G(O) 



fS(t-t>) 

(i-xfXi-x, 2 ) 

N 



5{t-t') 



(38) 
(39) 



while other covariances are of order a/r or smaller, see Appendix T5[ The jo int stochastic dynamics of allele frequencies 
and the correlation between loci has been studied by Ohta and Kimura ( |1969[ ) using a two-locus model. Here, we 
study a multi-locus model making the simplifying assumption that the recombination is faster than all other processes. 

In this case, the 2nd order cumulants relax much faster than allele frequencies change and we can solve the equation 
for xij assuming fixed X i- The solution can be decomposed into a deterministic component due to the competition 
between epistatic selection and recombination and a stochastic component. 



Xij(t) 



- Sxi 



(40) 



The deterministic component is the familiar QLE value from Eq. (27), while the stochastic component Sxtj has an 

(i-x?)(i-x?) 

± — / 



D 



We will now use this result to study 



auto-correlation (5xij(t)5xij(t + At)) — — A J,^ A ^ e - rAt ; see Appendix 
Langevin equation for X i- We have to distinguish the case where the deterministic component to Xij dominates over 
the stochastic term or vice versa. In order to compare the stochastic to the deterministic term, we have to average 
the former over the time scale of the dynamics of Xi given by the inverse of d ^ ~ /,. Recalling that the equilibrium 

value of the Xi ~ 1 — (t/fh we nn d that the deterministic contribution to Xij dominates if Nfi 3> 1 and fij 3> /i. In 
the opposite limit, the stochastic contribution 5xij will affect the dynamics of Xi more strongly than the deterministic 
one. We will now show how the equilibrium distribution of allele frequencies is affected by correlation between loci in 
these two cases. 
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A. Wright's equilibrium in the QLE approximation 



Assuming we can neglect the stochastic contribution to Xiji the Langevin equation for the Xi (interpreted in the 
Ito sense) corresponds to the following forward Kolmogorov equation for the dynamics of the probability distribution 
of allele frequencies by Q{{xi},t) (Gardiner 2004) 



d t Q({Xi},t) = J2 d x 



1 

2N 



E d Xi (XijQttXi}, *)) + Q({Xi}, t) 2nxi - E Xiid Xi (F) 



(41) 



This multilocus version of the diffusion equation for allele frequencies in linkage equilibrium (no correlations) appears 
already in Kimura ( 1955 1. It has a steady solution where all probability flux vanishes, i.e. where the term in brackets 



is zero for each i. In complete linkage equilibrium, the matrix Xij is diagonal and different allele frequencies decouple. 
One obtains the equilibrium distribution 



Q({xJ) = Ce 2 ^te»n(l-X?) 



2\2jVjU-l 



(42) 



where F({xi}) is the mean fitness evaluated in linkage equilibrium obtained by replacing each by its moment Xi in 
Eq. (Jlj). The term e 2NF ^ Xi ^ is analogous to the contribution of energy to a Gibbs measure, while Ili(l — X?) 2 ^^ 1 
plays the role of an entropy. Note that for 2N(i < 1, the distribution is singular at |x-i| = 1. In the opposite case 
2Nfi > 1, Q{{xi}) vanishes if any of the \xi\ = 1. Instead Q({xi}) has a maximum in the interior of the hypercube 
defined by \xi\ < 1- 

The corresponding solution for QLE, where Xij has small but steady off-diagonal entries is derived in Appendix [D| 
with the result: 



Qtixi}) = c, 



(43) 



The genotype distribution assumes this exponential (Boltzmann) form ~ e N( " F ^ since the mobility matrix Xij is pro- 



portional to the auto-correlation of the genetic drift. Eq. (43) provides a systematic extension of Wright's equilibrium 
to QLE, which appears to be a new result. 



B. Wright's equilibrium with stochastic linkage disequilibrium 



In absence of epistasis or in cases where selection is weak or comparable to the strength of genetic drift (diffusion 
constant), the deterministic expectation for Xij is small compared to its fluctuations. The coupling between different 
allele frequencies in Eq. ( 36 1 has therefore fluctuating sign and acts as an additional noise source with auto-correlation 
time (rcjj) -1 . Such an increased noise level increases the diffusion constant in the Fokker-Planck equation for each 
of the Xi by a factor 



N 



x 2 ) 



d(F)Y 



dx 3 J 



(44) 



This increase in diffusion constant is often phrased as a reduction in effective population size N e and is known as a 
manifestation of the Hill-Robertson effect (Hill and Robertson 1966). Note that the correction has the structure of 



the additive variance in fitness where each term is compared to the square of the recombination rate between the loci 



i and j. This result was derived in the context of fixation probabilities of novel mutations in Barton (1995b). It has 



been shown that this effective increase in the diffusion constant through stochastic correlations of loci can select for 
increased recombination rates (Barton and Otto 2005). 



C. Equilibration towards a steady state 



The approach to the equilibrium distribution is governed by the smallest non-zero eigenvalue of Eq. (41 1. For 
2N/J, > 1 and smooth fitness landscapes, this relaxation rate is governed by the larger of fi and the scale of selection 
on individual loci. The corresponding time scales can be very long. Furthermore, if different parts of sequence 



13 



space are separated by fitness valleys ("energy barriers"), relaxation to the steady state can take exponentially long 



(Weissman et al. 2010) 



Similar equations for the distribution of allele frequencies apply in the context of spatially structured populations, 
in which case the role of mutation is played by migration of individuals. The latter problem was the subject of work 



by Wright (1932). Migration rates and the associated influx of foreign alleles are often much larger than mutation 



rates and rapid equilibration is plausible. 



VII. BREAKDOWN OF QLE 



QLE greatly simplifies the dynamics of the genotype distribution, but the perturbation theory nature leaves one 
with the question about the range of its validity. In particular, we know from Statistical Physics that Gibbs measures 
of the form of Eq. (17) can lead to a so-called glass transition where the structure of the distribution changes 



qualitatively. Below the glass transition, different realizations of the system have a non-vanishing probability to be 
(largely) identical, which is quantified by the overlap distribution (Parisi order parameter, |Mezard and Montanari 
(2009)). A related transition in which the population condenses into a small number of genotypes is driven by the 
competition between epistasis and recombination. It occurs already in the deterministic mean field setting and is 
discussed below in Section VILA QLE can also become unstable at low recombination rates even in the absence of 



epistasis because of the discreteness of contributions of individual loci in a finite genome. The instability in that case 
is driven by fluctuations due to finite population size and is discussed in Section |VII.B| 



A. Infinite TV and L limit: Alleles vs. Genotypes 

To gain some heuristic insight into the range of validity of the perturbation expansion in cr/V, it is useful to 
study the following coarse-grained "quantitative genetic" version of QLE which yields an explicit criterion for the 
validity of QLE (iNeher and Shraiman 2009). Instead of following the entire genotype distribution, consider the joint 
distribution P(A, E, t) of additive A and epistatic E contributions to fitness defined via A = Fa{q) (comp. Eq. (30)) 
and E = F — A. Hence additive and epistatic contributions are defined with reference to the current distribution of 
genotypes. The joint distribution of A and E evolves according to 



d t P{A, E,t) = {A + E- (A) - (E))P(A, E,t)+r (p(E) J dE 1 P(A, E', t) - P(A, E, t)\ 



(45) 



where (A) and (E) are the mean additive and epistatic fitness in the population. Here we have assumed that the 
epistatic fitness of novel recombinants is independent of its parents and given by a random sample from the density 
of possible epistatic fitness values p(E) (the "house-of-cards" model, ( Kingman] 1978)). We assume p(E) to be a 
Gaussian with the variance equal to <rj - the epistatic component of fitness variance defined in Eq. (31). Additive 
fitness of recombinants is a random sample from the current distribution of additive fitness in the population, i.e. the 
marginal J dE' P(A, E' , t). The model does not include finite population size effects and assumes that both A and E 
are from a continuous distribution. The latter implies that the number of loci L that contribute to fitness is very large, 
while the individual contributions of loci are small (comp. Section [H]) . In this sense, it is a deterministic mean-field 
model. 



Equation (45) has a factorized solution P(A,E,t) = 9{A,t)u)(E) with 



6{A,t) 



(A- (A))* 



where 



dt 



(A) 



(46) 



u(E) 



rp(E) 
(E)-E 



where (E) = J dE Euj(E) 



The mean epistatic fitness is determined by enforcing the normalization of u)(E), i.e. J dEui(E) — 1. Note that this 
solution is a QLE solution: Fitness increases with a rate given by the additive variance, while the epistatic contribution 
to fitness is steady with a magnitude controlled by recombination. Unlike Eq. ( p7| , Eq. ( [46] ) implies a condition on r 
and the density of states: p{E) has to vanish for E >r+ (E). Otherwise, oj is not normalizable. The density of states 
p(E) is typically of Gaussian form, and given 2 L states has its maximum at E max w <rjy/2L log 2 (if N <C 2 L , as will 
be generically the case, E m , 



(J]\/2 logN). Hence QLE is expected to break down at r c « E n 



(E) « E n 
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FIG. 2 Genetic interactions and the breakdown of QLE. Panel A shows the range of validity of QLE as a function of a/r and the 
heritability, i.e. the ratio of additive variance to the total fitness variance. Below the transition line, strong linkage equilibrium 
is expected and selection operates on genotypes rather than alleles. Panel B shows a illustration of a possible "interaction 
graph" of polymorphisms on a block of chromosome. Loci interact with nearby loci, as well as with distant loci outside the 
block. The epistatic fitness variance solely within the block, i.e. averaged over the rest of chromosome, is proportional to the 
number of interaction terms (arcs in the figure). 



The dynamics of the distribution of P(A, E) change dramatically as r falls below r c : For r > r c no genotypes with 
E > r + (E) exist. Hence all genotypes are destroyed by recombination and short-lived. At r < r c , however, many 
genotypes with E > r + (E) exist which can outrun recombination and grow exponentially. The genotype distribution 
is no longer a product of additive and epistatic parts, but contains clones which are populated by many individuals. 
Selection now operates on the entire genotype over many generations and the relevant dynamical quantities are now 
clone sizes rather than allele frequencies, which are slaved to the performance of the clones. The alleles that make up 
the most successful genotype will fixate, not necessarily those with the most favorable additive effect. The transition 
line between the two regimes is sketched in Fig. [2] with the ratio of recombination to selection on the x-axis and the 
"heritability" - the ratio of the additive variance to the total variance h 2 — o\/{cr\ + erf) - on the y-axis. (Note 
that heritability also measures the correlation between fitness of a recombinant offspring and parental mean |Lynch| 
and Walsh (1998).) At low recombination rates and strong epistatic interactions, selection operates on genotypes 



while at high recombination rates or in absence of epistasis, selection operates on the additive effects of alleles. The 
distinction between genotype and allele selection regimes goes back to Franklin and Lewontinj ( 1970 ); Slatkin (19721, 
who showed that a related transition occurs in models with strong heterozygote advantage. The regimes of allele and 
genotype selection are summarized in Fig. [2] In absence of epistatic interactions or heterozygote advantage, a similar 
condensation phenomenon occurs only at very low outcrossing rates r = 0(N~ l ) ( jRouzine and Coffin[ |2005[ ) . 

The condensation of genotypes goes along with a dramatic speed up of the allele frequency dynamics: In QLE (allele 
selection) each allele frequency is driven by an effective additive coefficient ai ~ a/V~L (each ai accounts for ~ L _1 
of the additive variance a 2 A ). When selection operates on genotypes, the time scale of selection is driven by fitness 
differences between individuals, which are of order a. The rate of chang e of allele frequencies is therefore greater by 
a factor \fh which could be a large effect (Neher and Shraiman 20091. We emphasize that the stationarity of the 



distribution of the epistatic component of fitness u>(E) holds only on time scales short compared to that of the allele 
frequency dynamics. 

The simple picture of the transition in a facultatively outcrossing species can also apply to blocks of chromosome in 
obligate sexuals. Consider a block that harbors / loci spread over a map distance c (on average c recombination events 
within the block per generation). If epistatic fitness within the block exceeds c, QLE will break down since individual 
haplotypes will be amplified by selection above less fit recombinants. Whether such local breakdown of QLE will 
occur depends on how fitness variance and recombination rate depend on the block size. The recombination rate is 
proportional to the block size, and, assuming constant density of polymorphism, will be proportional to the number I 
of polymorphic loci. Similarly, the additive variance is proportional to I and the r-m-s therefore ~ yl. The epistatic 
variance within the block scales with the number of interactions between loci within the block, as illustrated Fig. [2}3. 
Any given locus will interact only with a fraction of all other loci, i.e. fij is sparse, and the number of interactions 
between loci within a block depends on whether these sparse interactions tend to be local or not. If any two loci 
are equally likely to interact, the number of interactions within the block is <~ l 2 , so that r-m-s epistatic fitness is 
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~ I ~ c. Hence the ratio of recombination within the block and the epistatic fitness are independent of the block 
size and QLE is either globally stable or unstable. A different conclusion is reached if interactions are local and each 
locus interacts with k nearby other loci. As before additive fitness ~ y/l, but the number of interactions within the 
block is ~ Ik. Hence the typical epistatic fitness is ~ \lk, which decreases less fast than c as the block length is 
decreased. We therefore expect that QLE is unstable on scales below a critical block size l c , where local epistasis 
overwhelms rare recombination. This local selection on coadapted haplotypes can coexist with establishment of QLE 
on longer genomic scales (Neher and Shraiman 2009). We will come back to the recombination and selection on 
different chromosomal scales in the Discussion. 



B. Validity of QLE for finite N and L. 



The above discussion of the breakdown of QLE focussed on the competition between genetic interactions driving 
and recombination destroying correlations in the limit where fluctuations are negligible and contributions of individual 
loci are small. We will now discuss how the discrete contributions of individual loci and the number fluctuations in 
finite populations can drive populations off the QLE manifold. This problem has a long history in population genetics 
and was mainly discussed for scenarios without genetic interactions, i.e. on the line where the heritability equals 1 in 
Fig. [2pV. In this case the only source of correlations are the initial condition or fluctuations. Maynard Smith ( 1968 ) 
showed that without genetic interactions and with no correlations in the initial condition, correlations do not develop 
in an infinite population at any recombination rate, in accordance with Fig. [2]A However, novel mutations arise in 
single copies on random genomes, giving rise to correlations: QLE has to be stable with respect to these perturbations. 

The hallmark of the QLE approximation are slowly changing allele frequencies and steady and perturbative cor- 
relations between loci. The latter will only be true, if the correlations relax, i.e. are governed by an equation of the 
form Xij = P — a Xij with a — 2(fiXi + fjXj) + rc ij > (ignoring mutations). Hence the QLE state is unstable if 
—2{fiXi + fjXj) > rc ij- In that case any small deviation from Xij = 0, which could be due to stochastic fluctuations, 
will grow. This effect has important implications for the evolution of recombination: Consider two closely linked loci 
at which beneficial mutations happen. Both novel mutation exists initially as a single copy (xi ~ —1) and will most 
likely reside in different individuals, i.e. are anti-correlated or in negative equilibrium. Selection will now amplify the 
initial Xij if 2/, + 2fj > rc^, generating predominantly negative LD. This growth of correlations due to selection on 
individual loci slows down adaptation and can result in the loss of beneficial alleles. This phenomenon is known as 
Hill- Robertson interference and it contributes to potential benefits of sexual reproduction (Barton 1995b Barton and 



Otto 2005 Hill and Robertson 19661. While this cumulant based approach to interference between sweeping loci is 
tractable for few loci, it becomes intractable in populations in which many sweeping loci are tightly linked (Cohen 
et al 



2005 Neher et al. 2010 Rouzine and Coffin 2005). 



C. Cumulant analysis beyond QLE 

Even though the QLE approximation breaks down when correlations are no longer slaved variables, the cumulant 
expansion can be useful to study the short term dynamics of systems with a small number of loci, in particular if the 
initial conditions are such that higher order cumulants are small. Furthermore, if only a few isolated pairs of tightly 
linked loci are present, cumulants between these pairs can be treated as dynamical variables, while all other pairs 



for which the Xij are stable are treated in QLE. Such an analysis has for example been performed by Stephan et al. 



(2006) to study LD patterns between neutral markers following a selective sweep. 

Explicit modeling of stochastic multi-locus systems typically requires computer simulations, which are computa- 
tionally expensive when the number of loci or the population size is large. However, making use of the Fast-Fourier 
Transformation on the 2 L dimensional genotype space, one can speed up such simulation from a runtime that scales 
as 8 L to 3 L . The FFT allows to calculate and reuse the frequency of subsets of loci from which the distribution of 
recombinant genomes can be assembled. An efficient implementation of multi-locus evolution for arbitrary fitness 
functions and genetic maps is available from the author's website. Cumulant equations to higher order involve "book- 
keeping" of many terms and is best done with computer algebra systems. A package for Mathematica® has been 
developed by |Kirkpatrick et aL| ( |2002| . A implementation for Maple® is available from the authors. 



VIII. DISCUSSION 



We have presented a review of the dynamics of multi-locus genotype distributions and the resulting dynamics of 
quantitative traits. We focused in particular on how the distribution of genotypes can be parameterized by allele fre- 
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quencies in the weak selection/fast recombination limit. This description extends beyond "beanbag genetics" allowing 
also for weak correlation (i.e. linkage disequilibrium) between loci. The central element is the Quasi-Linkage Equi- 
librium approximation pioneered by Kimura. QLE emerges as a perturbation expansion in the weak selection/rapid 
recombination limit similar to high-temperature expansion in statistical physics. In a suitably defined system, the 
population genetics can be classified by the ratio of the strength of selection and the rate of recombination and the 
degree to which the fitness variation is additive or epistatic, see Fig. [2] At high recombination and additivity, QLE is 
an accurate approximation. This regime is separated from a regime at low recombination and strong epistasis, where 
QLE breaks down and the population condenses into a few fit genotypes. 

Our exposition assumes a panmictic, random mating, and haploid population. While the former are pretty common 
assumptions, assuming haploidy in recombining population might raise objections. Our aim was to discuss the 
interplay between selection, genetic interactions and recombination in multi-locus systems. Dominance is a special 
kind of genetic interaction, where a locus interacts with itself, giving rise to additional non-linearities. These non- 
linearities can stabilize loci at intermediate allele frequencies, a process not possible in haploid populations. The effects 
of dominance, however, are well understood at the single locus level, as well as when many loci with heterozygote 



advantage are close to each other (Franklin and Lewontin 1970). Within QLE, the dynamics of allele frequencies in 



diploid populations is still relaxational and maximizes the mean diploid fitness. A full parametrization of the diploid 
populations and diploid fitness requires a straightforward, if somewhat tedious, generalization: To represent diploid 
one should (i) double the number of loci, (ii) define a genetic "transfer function", C({£}), that represents meiotic 
crossover of the parental genomes, (iii) extend the fitness function F(g) to 2L hypercube to parameterize the 3 L states 
(homo- and heterozygocity at L loci). Another simplification of our exposition was the use of the continuous time 
description, in contrast to the more common discrete generation formulations of population genetics. Continuous time 
formulations assume that the population changes little in one generation. If this is the case, the results are completely 
equivalent and one can make use of calculus instead of recursions and difference equations. 

The QLE approximation will often be appropriate for panmictic populations where genetic variation is replenished 
by de novo mutations. In this scenario, novel mutations establish if they blend in well with genetic make up of the 



population. This is manifest in the QLE equation (Eq. (27)) were alleles are selected on the basis of their additive 



effect, i.e. their effect on fitness marginalized over the distribution alleles at other loci in the population. The fitness 
of individual genotypes is not relevant to the evolutionary dynamics, since genotype frequencies are determined by 
allele frequencies (and sampling noise in finite populations). This issue was recently discussed in Livnat et aL|(|2008|). 



A very different evolutionary dynamics follows a hybridization event (Barton 2001| Nolte and Tautz^ 2010[ |Orr[ 
1995[ ), i.e. a situation when two strains of one species that have been evolving in isolation for some time come in contact 
again. The two strains will differ at many loci and these differences have never been tested for compatibility. Crossing 
two such diverged strains can result in a phenotypically diverse population from which novel hybrid species can emerge 
( Nolte and Tautz||2010[ ). Such speciation after hybridization is similar to the clonal population structure observed in 



theoretical models of the selection dynamics after hybridization (Neher and Shraiman, 2009). In this regime of clonal 
competition, the allele frequencies are slaved to the dynamics of the clones and the average effect of an individual 
mutation affects the fate of a clone only very mildly. The lucky accident that produced through recombination a 
very fit genotype that contains the allele determines whether the allele can fixate or not. The possibility of a sharp 
transition between mixing and not-mixing of two populations in a hybrid zone has already been described by |Barton 



(1983). who used a model of hybrid-inferiority. In the limit of large number of contributing loci, there exists a critical 
ratio of recombination rate and selection against hybrids, which separates the regimes of mixing and non-mixing. 
Similar transitions are expected if the reason for out-breeding depression is epistasis rather than dominance. 

The qualitative differences between the genotype and allele selection regimes also sheds light on the importance of 
stochasticity (genetic drift). Allele frequencies are well sampled by O(N) copies, unless the allele is very young (or 
about to go extinct). Stochasticity therefore matters only during the establishment phase of the allele. As soon as the 
frequency exceeds (AT/j) -1 ps y/L/Na, selection dominates. In the genotype selection phase, however, the founding 
of each genotype can, if it is exceptionally fit, change the fate of the population dramatically. 

The transition to genotype selection driven by epistatic interaction is related to spin-glass transition in models for 
disordered physical systems and magnets. Within these models, the probability of finding the system in a particular 
state {s^ is given by 

P({si}) ~ e - n (i s i})/ kT = e -w[l>i*<+£« ^ ( 47 ) 



and hence completely analogous to Eq. (17). Such system generically reside in one of three states: paramagnetic, 
ferromagnetic, glassy. In the paramagnetic state at high temperature different parts of the system are uncorrelated, 
which is analogous to QLE. The perturbation expansion in a/r is very similar to a high temperature expansion in 
statistical physics. At low temperature, the behavior depends on the structure of the Hamiltonian H({si}). If most 
of the Jij have the same sign, the system will go to an energetically favored ordered state where spins are aligned, 
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giving rise to a ferromagnet. In this case "H({si}) has one heavily preferred energy minimum, corresponding to a very 
fit genotype. 

A different low temperature behavior is found when the have erratic sign. In that case, not all interactions 
can be in their favorable state simultaneously and the resulting landscape has many minima and maxima. At low 
temperature, the system condenses into one of the minima. This "spin-glass" phase is characterized by a non-trivial 
overlap distribution: Different realizations of the system, drawn from the ensemble defined by Eq. (47), will fall into 
clusters of different degrees of similarity (measured by Hamming Distance). The clusters themselves have subclusters, 

1987[ ). This is in contrast to the high temperature 



giving rise to a hierarchical ultrametric structure (Mezard et al. 



phase, where most systems share a typical number of sites. These qualitatively different overlap distributions above 
and below the spin-glass transition have a direct analogies to population structure and heterozygosity: In the high 
recombination limit, genotypes in the population are assembled from the available alleles more or less at random 
such that any two individuals differ at about 2j^f , z^(l — Vi) sites (z^ being the allele frequency at locus %). At low 
recombination (or substantial inbreeding), the population will condense into fit genotypes (or inbred groups) that are 
much more similar to each other than to members of the general population. 

Figure [2^3 illustrates pair- wise interaction between polymorphic loci along the chromosome. In general, we expect a 
complex and possibly hierarchical pattern of interactions: A given pair of distant genes will have only a low probability 
to interact substantially, while polymorphisms within one gene and its regulatory elements are much more likely to 
interact strongly. Nearby polymorphisms in a protein (Callahan et al. 2011) will be still more likely to interact. In 



obligate sexuals, the sparse long range interactions will rarely suffice to produce appreciable correlations between loci. 
Within small stretches of chromosomes, however, recombination rates are low and if the strength of interactions within 
this stretch is sufficiently high, QLE will locally break down. Consider for example a one centimorgan long region, 
which in humans corresponds to about one mega base and harbors around a thousand polymorphisms. If the typical 
epistatic contribution fitness of this stretch of chromosome were on the order of 1%, we expect run-away selection on 
coadapted haplotypes and strong correlations. Since distant parts of the genome are in QLE, one expects a "module" 
selection regime, where loosely linked and weakly interacting "modules" are in QLE, but strong interactions and 
infrequent recombination has led locally to condensation into coadapted haplotypes ( Neher and Shraiman| 200 
(An excellent early discussion of such epistasis driven "coagulation" in the "soup" of genes is found in (ITurner 



1967).) Put otherwise, we can view such a system as consisting of weakly-interacting mesoscopic loci, at which several 
These super-alleles are "destructible" , in the sense that recombination within leads to reduced 
However, since recombination within these alleles is rare, quantitative traits would 



super-alleles segregate, 
fitness and purging by selection, 
be highly heritable on short timescales and quantitative genetics would work as usual. 

Our discussion of the multilocus theory and QLE was guided by ideas of Statistical Physics. The explicit form of 
the (approximate) genotype distribution function P(g) parametrized by instantaneous allele frequencies is the central 
pillar connecting the dynamics of population average traits - the subject of QG - to the individual-based evolutionary 
process. It is essential that the QLE distribution is reached on a relatively fast time scale of mating and recombination. 
Allele frequencies are well defined and vary slowly on this time scale. The QLE ensemble should not be confused with 
a very different mutation/selection/drift ensemble - which could be rightfully termed the "Wright equilibrium" (Eq. 
(42 1) - which is often invoked as a link between evolutionary dynamics and Statistical Physics. Wright equilibrium gives 
a stationary distribution of allele frequencies which would be established in a finite population (N playing the role of 
inverse temperature) on a time scale longer than the inverse mutation rate fx, provided stationary selection pressures. 



Ruggedness of the fitness landscape could further increase this equilibration time scale exponentially ( jWeissman et al. 



2010[ ). Clearly, this type of equilibrium applies on a very different time scale than the phenomena addressed in the 



present work. Related ideas were developed in the context of quasi-species theory to study the conditions under which 
hereditary information can be maintained over long times ( |Eigen 1971 Franz and Peliti 1997). The focus of these 
studies was pre-biotic evolution, where fidelity of replication was most likely low and stability of genomic information 
can be sens ibly studied using a simple equilibrium model. Equilibrium argume nts were also applied to evol ution of 
codon bias (Iwasa 1988) and the evolution of transcription factor binding sites (Mustonen and Lassig 2005). In the 



latter two cases, an ensemble can be constructed by combining many instances of the same sequence motive which 
was under constant selection pressure for very long time (conserved transcription factor binding motive or conserved 
preference of certain codons over others). In many cases, however, the equilibrium state is of little relevance. 

In conclusion, in this review we have provided a derivation of the genotype distribution in the QLE approximation, 
providing a systematic generalization of Fisher's theorem, Kimura's diffusion theory and Wright's equilibrium from 
LE to QLE, which includes the effect of (weak) correlations between loci. We have also discussed the limitation of 
the QLE approximation and the structure of the genotype distribution at low recombination rates. 

It is our hope that better understanding of the QLE approximation will promote progress in understanding the 
effects associated with its breakdown, whether due to strong epistasis or strong physical linkage, such as for example 
the Hill-Roberson effects (hitch-hiking and background selection) which still await comprehensive treatment. 
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Appendix A: Glossary 

Allele State of a locus, for example the base A, C, G or T at a certain position 

Crossover rate In meiosis, parental chromosomes are paired up and crossed over. The density of crossovers on the 
chromosome is called crossover rate. 

Dominance Interaction of the two alleles at the same locus in diploid organisms 

Epistasis Genetic interactions between alleles at different loci, i.e. a dependence of the effect of an allele at one locus 
on the remainder of the genome. 

Fitness Expected reproductive success of an organism. For modeling purposes, this is often equated with the growth 
rate (Malthusian or log fitness) or the average number of offspring in the subsequent generation (absolute fitness) . 

Gametes Egg and sperm 

Genetic Drift Sampling fluctuations of genotype or allele frequencies. Genetic drift enters as the diffusion term in 
the Fokker-Planck equation for the dynamics of the distribution of allele frequencies. 

Genetic map The cumulative crossover rate along the chromosome. The average number of crossover events per 
chromosome is the map length. 

Genotype State of the genome, i.e. the set of alleles an individual carries. 

Haplotype Alleles inherited from one parent. In diploids, two haplotypes make one genotype. 

Heritability Broad sense heritability is the genetic component of traits, i.e. the concordance of traits between 
monozygotic twins. Narrow sense heritability refers to the genetic component of traits that is inherited in 
sexual reproduction, i.e. the correlation between trait values of parents and children. 

Heterozygosity Fraction individuals in a diploid population that carry distinct alleles at a locus. 

Homozygosity The complement of heterozygosity 

Linkage Loci on the same chromosome are linked and share history until crossover events separate them. 
Linkage (dis)equilibrium Absence (presence) of correlations between loci 
Locus Location on the chromosome, e.g. a gene. 

Mean Fitness To preserve overall population size, fitness is often measured with respect to the mean fitness of the 
population. 

Meiosis Division of a diploid cell to produce haploid gametes 

Panmictic A population is panmictic if each individual is equally likely to compete and interact with any other 
individual. In practice, this requires that dispersal is fast compared to population genetic time scales. 

Polymorphism A locus with variation, i.e. the population contains several alleles at this locus. 

Random mating Simplifying assumption that mating is independent of genotype, phenotype, and environment. 

Recombination Process of reshuffling of the genetic material in sexual reproduction. 

Outcrossing Fertilization with sperm/pollen from a different individual 

Selfing Many plants and other organisms have female and male sexual organs and can self-fertilize or self-pollinate. 
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Appendix B: Notation 

The specification of genotypes and parameterization of genotype-phenotype maps is not unique and our notation 
differs from the traditional population genetics choice. Conventionally, one chooses one "wild-type" reference genome 
(0,0,..., 0) and enumerates deviations from this reference. This is useful when a well defined wild-type genotype exists. 
In diverse populations, for example the progeny of cross between diverged strains, the reference free parameterization 
we are using here is more natural. The allelic state at each locus is denoted symmetrically by ±1, e.g. whether an 
allele comes from one or the other strain. The two different parameterizations are completely equivalent and related 
to each other by a simple linear transformation (see Table |H| below). In the present context the reference free notation 
simplifies the algebra since the Sj = ±1 basis is orthogonal when averaging over the genotype space. The relation to 
the Fourier transform allows an unambiguous decomposition of the fitness function into additive parts and epistatic 
components of different order (Parceval's Theorem), while in the reference based parametrization of fitness functions, 
more akin to a Taylor expansion, coefficients depend explicitly on the choice of reference. We have also deviated from 
the traditional notation for linkage disequilibrium because we want to use the diagonal \u — 1 — X? components 
of the cumulant matrix (two times the heterozygosity at locus i) on the same footing as the off-diagonal one. 



Symbol 


Meaning /Definit ion 


9 


Haploid genotype: g = {si, . . . , sl} 


P(g,t) 


Genotype distribution in the population 


(...) 


Population average 


F(g) 


Fitness (growth rate) of genotype g 




Contribution to fitness of the i\ . . . ik set of loci 




Additive effect of locus i 


„2 2 2 

a , a A , o i 


Total, additive, and epistatic variance in fitness 


ii G {0,1} 


Origin of locus i, i.e. maternal or paternal 


C({£}) 


Probability of the recombination pattern 


Cij 


Probability that loci i and j derive 




from different parents 


/i, r 


Mutation and outcrossing rate 



TABLE I Table of symbols 



Quantity 


Our notation 


Allele at locus i, {a,i,Ai} 




Si G {-1,1} 


Allele frequency Vi 




v% = (1 + xO/2 where %i = ( s >) 


Linkage Disequilibrium Dij 


(i ^ j) 


4Aj = Xij = ( s i s j) ~ ( s i)( s ]) 


Heterozygocity Hi = 2^(1 - 


- Vi) 


2Hi = xu = (1 - Xl) 



TABLE II Population genetic quantities in our notation 



Appendix C: QLE in terms of effective fields 

In this appendix, we discuss how the fields <f>i and <pij introduced to parameterize the genotype distribution P(g, t) 
in Eq. (17) are related to the cumulants of P(g,t). We also detail how the recombination term in Eq. Q can be 



evaluated explicitly within the QLE perturbation theory. We parameterized the genotype distribution via 

log P(g,t) - $(t) + Y,<l>i(t)si + J2 ( t ) v( t ') SiS j- ( C1 ) 



i<j 



The constant term is determined by the normalization of the distribution, the coefficients <pi (t) are related to frequen- 
cies and the second order coefficients 4>ij(t) to the connected correlation between loci. In the limit under consideration, 
the second order contributions are small and we evaluate the coefficients to leading order in <f>ij{t). 
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9 



(C2) 



2 L 1 + E tanh(0fc) tanh(</>j) cosh (</>,;) 
fe<j 



The relations between x%i Xij an d 0i ; 4>ij given in Eq. (20) follow by differentiation. 

To arrive at the equations for the time evolution ofthe fields fa and <f>ij (Eq. ( 24 1 ) , we have to evaluate the 
recombination term in Eq. (23 1. This is done below. The terms proportional to fa(t) cancel exactly between numerator 
and denominator and we arc left with 



E wftW) 

{£»}{<} 



P( g W)P( g (/)) 

P{g)P{g') 



r E <7({f}W) 
{&}«} 



e Ei<j ^«[K<Si+f<&-i)(»<*j+*'i4)+(^&+f«^)( a ' s j+ a i^)] - 1 



(C3) 



In the limit under consideration, the second order contributions have to be small enough that the entire exponent is 
small. In this case, the exponential can be expanded and the different terms averaged individually. 



E Cm)P(g') 
{€>}«} 



P{g (m) )P{g U) ) 



- 1 



P{g)P{g') 

E E hi + Siti - WW + + + ^')(«i(*j> + ( s i) s j)] (C4) 



where c« is the probability that an odd number of crossovers happened between loci i and j. 



Appendix D: Diffusion theory and Wright's equilibrium 



In this appendix, we detail intermediate steps to arrive at the diffusion equation for the allele frequencies in QLE 
and the generalized Wright equilibrium. The noise terms in the Langevin equation ( 36 1 stem from the multinomial 
sampling of the genotypes or gametes. From the covariance of the multinomial distribution, we can therefore determine 
the covariance of the noise terms Q for the Xi an( i the for the Xij- The covariance the changes in Xi an d Xi> f° r 
example, can be calculated as follows 



gg' 



1 

N 



E s iSj p( g )(i - P(g)) - E ^P(g)P(g') 

a g+g' 



N 



(Dl) 



The other covariance terms can be calculated analogously. In particular, one finds (A X 

x?)(i-x, 2 )- 



N~ 



iX jj =N-'(l 
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The effect of deterministic correlations 



If deterministic correlations dominate over the fluctuations in Xij the forward Kolmogorov equation for the distri- 
bution of the Xi is given by 



dtQ({Xi},t) =J2d x 



^ E d xi (xijQdxi}, *)) + Q({xi}, t) \2fiXi - E < F > j 



(D2) 



In the steady state, all probability fluxes vanish. The i component of the probability flux is precisely the expression 
in brackets above and hence has to be equal to zero. Multiplying the bracket with 2N\~^l and summing over i (XT. 
is the ki element of the matrix inverse of Xij)> we have 



d k Q({xi}) = Qdxi}) - E x^djxa - £ XuXi + d Xk (f) 



(D3) 



Next, we use the fact that the off-diagonal elements of Xij are small and Xij — 7y(l — Xi )(1 — Xj)> while Xu = 1 — Xi ■ 
To first order in the off-diagonal elements, the inverse is given by x]±_ = (1 — X?) -1 an d off-diagonal elements 



X^ = — Xij(l — xf) (1 — Xi) — ~7ij- Going over the terms in Eq. (D3| one by one, we have 
E *ki d XiXij = E *W E d xsXij + E Xkl d XiXii 

ij i j^ti i 

= E XuXij E " x]) + E XuXiidxt log(l - X?) 

i jjti i 

= d x k M 1 - Xk) 
The mutation term can be evaluated as follows 

m v E XklXi = ~2N»d Xh log(l - x') - 4^ E T«Xi 



Substituting these terms into Eq. ( D3 ) and 7^- 



-2Nfid Xk log(l - Xk) - ±Nfid Xh JikXiXk 

i^k 

— , we have 



(D4) 



(D5) 



d k Q({Xi}) = Q({Xi})d Xk I (2iV M - 1) log(l - xl) + 27V[(F) + 2^ ]T 



i^k 



rc ik 



which is straight-fowardly integrated to 



Q({ Xi }) = Ce 



2iV<F)+4JV A i£ i<; 



fikXjXk 



W-xl 



2\2Nfj,-l 



i=l 



(D6) 



(D7) 



The effect of fluctuating correlations between loci 

Even when associations between loci are zero on average, fluctuations of Xij can affect the allele frequency dynamics. 
The coupling between different loci acts as an additional noise source on the dynamics of allele frequencies. Grouping 
deterministic and stochastic forces, the corresponding Langevin equation for Xi is given by 



Xi (t + At) - Xi(t) = At [ Xll d Xl (F) - 2/%] 



t+At 



Ex«(t')^.<F> 



(D8) 
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where the integral constitutes the fluctuating noise term. Solving the Langevin equation for \ij (t) assuming constant 
Xi and Xj: one finds 



(1-X?)(l-X?)e-^ Z 



2Nrc. 



(D9) 



Averaging the square of the noise term in Eq. ( D8 ) , we find 



t+At rt+At 

dt' I dt" 
t Jt 



J2xikd Xj {F)+Q 

k^i 



X»AT 
N 



(D10) 



The cross-term is of order N 3 ^ 2 and can be neglected. 
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